MODULE sefmm_driver

   USE fmm_global_paras
   IMPLICIT NONE
   PRIVATE
   ! Public procedures
   PUBLIC :: sefmm_get_GTO_FF_J_matrix,    &
             sefmm_get_FF_J_matrix,        &
             sefmm_get_NF_J_matrix

CONTAINS

!-------------------------------------------------------------------------------

   SUBROUTINE sefmm_get_GTO_FF_J_matrix(scheme,basis,dens,J_matrix)

      USE fefmm_box_utils,          ONLY: fefmm_init_box_utils,           &
                                          fefmm_get_boxes
      USE fmm_shell_pairs,          ONLY: fmm_get_shell_pairs
      USE fmm_driver,               ONLY: fmm_get_multipole_potential
      USE fmm_boxed_multipole_ints, ONLY: fmm_pack_boxed_mpoles,       &
                                          fmm_contract_multipole_potential

      IMPLICIT NONE
      TYPE(scheme_paras), INTENT(IN)    :: scheme
      TYPE(fmm_basis),    INTENT(IN)    :: basis
      REAL(REALK),        INTENT(IN)    :: dens(basis%nBas,basis%nBas)
      REAL(REALK),        INTENT(INOUT) :: J_matrix(basis%nBas,basis%nBas)

      TYPE(fmm_sh_pairs), POINTER :: sh_pairs(:)
      REAL(REALK), ALLOCATABLE :: Vpot(:,:)
      REAL(REALK)   :: dummy_dens(1,1) = one
      REAL(REALK)   :: T0, T1
      INTEGER(INTK) :: lmax, lmdim, nbox

      CALL CPU_TIME(T0)
      CALL fmm_get_shell_pairs(basis,sh_pairs)
      CALL fefmm_init_box_utils(scheme%grain,sh_pairs)
      CALL fefmm_get_boxes(nbox)
      IF (nbox < 3) RETURN

      lmax = scheme%raw_LMAX
      ! Pack moments in same box to reduce memory in FMM readin
      CALL fmm_pack_boxed_mpoles(lmax,basis%nBas,dens,nbox)

      ! Get multipole potential for each box
      lmdim = (1 + lmax)**2
      ALLOCATE (Vpot(lmdim,nbox))
      CALL fmm_get_multipole_potential(FE_FMM,dummy_dens,Vpot)

      ! Contract potential with AO moments box by box
      CALL fmm_contract_multipole_potential(basis%nBas,nbox,lmax,Vpot,J_matrix)
      CALL CPU_TIME(T1)
      CALL TIMTXT('>>> TIME USED get FF J-matrix', T1-T0, LUPRI)

      DEALLOCATE (Vpot)

   END SUBROUTINE sefmm_get_GTO_FF_J_matrix

!-------------------------------------------------------------------------------

   SUBROUTINE sefmm_get_FF_J_matrix(scheme,basis,Jmat,energy)

      USE fefmm_box_utils,       ONLY: fefmm_init_box_utils,           &
                                       fefmm_get_boxes
      USE fmm_shell_pairs,       ONLY: fmm_get_shell_pairs
      USE fmm_driver,            ONLY: fmm_get_multipole_potential
      USE sefmm_cheb_multipoles, ONLY: sefmm_build_cheb_multipoles,      &
                                       sefmm_contract_cheb_multipoles

      IMPLICIT NONE
      TYPE(scheme_paras), INTENT(IN)    :: scheme
      TYPE(fmm_basis),    INTENT(IN)    :: basis
      REAL(REALK),        INTENT(INOUT) :: Jmat(basis%nBas,basis%nBas)
      REAL(REALK),        INTENT(OUT)   :: energy

      TYPE(fmm_sh_pairs), POINTER :: sh_pairs(:)
      REAL(REALK), ALLOCATABLE :: Vpot(:,:)
      REAL(REALK)   :: dummy_dens(1,1) = one
      REAL(REALK)   :: T0, T1
      INTEGER(INTK) :: lmax, lmdim, nbox

      CALL CPU_TIME(T0)
      CALL fmm_get_shell_pairs(basis,sh_pairs)
      CALL fefmm_init_box_utils(scheme%grain,sh_pairs)
      CALL fefmm_get_boxes(nbox)
      IF (nbox < 3) RETURN

      ! Build multipole expansions of density for each box
      lmax = scheme%raw_LMAX
      CALL sefmm_build_cheb_multipoles(scheme,lmax,nbox,nbox)

      ! Get multipole potential for each box
      lmdim = (1 + lmax)**2
      ALLOCATE (Vpot(lmdim,nbox))
      CALL fmm_get_multipole_potential(FE_FMM,dummy_dens,Vpot)

      ! Contract potential with multipoles box by box
      CALL sefmm_contract_cheb_multipoles(nbox,lmax,Vpot,energy)

      CALL CPU_TIME(T1)
      CALL TIMTXT('>>> TIME USED get FF energy', T1-T0, LUPRI)

      DEALLOCATE (Vpot)

   END SUBROUTINE sefmm_get_FF_J_matrix

!-------------------------------------------------------------------------------

   SUBROUTINE sefmm_get_NF_J_matrix(scheme,basis,Jmat,energy)

      USE fmm_shell_pairs, ONLY: fmm_get_shell_pairs
      USE fefmm_box_utils, ONLY: fefmm_init_box_utils,           &
                                 fefmm_get_boxes
      USE sefmm_V_driver,  ONLY: sefmm_init_potential,           &
                                 sefmm_free_potential,           &
                                 sefmm_build_NF_potential
      USE sefmm_J_driver,  ONLY: sefmm_get_J_matrix,             &
                                 sefmm_contract_NF_potential

      IMPLICIT NONE
      TYPE(scheme_paras), INTENT(IN)    :: scheme
      TYPE(fmm_basis),    INTENT(IN)    :: basis
      REAL(REALK),        INTENT(INOUT) :: Jmat(basis%nbas,basis%nbas)
      REAL(REALK),        INTENT(OUT)   :: energy

      TYPE(fmm_sh_pairs), POINTER :: sh_pairs(:)
      INTEGER(INTK) :: nbox
      REAL(REALK)   :: T0, T1

      CALL CPU_TIME(T0)

      CALL fmm_get_shell_pairs(basis,sh_pairs)
      CALL fefmm_init_box_utils(scheme%grain,sh_pairs)
      CALL fefmm_get_boxes(nbox)
      CALL sefmm_init_potential(nbox,scheme)

      CALL sefmm_build_NF_potential(nbox,scheme%grain,scheme%NN_order)
      CALL sefmm_contract_NF_potential(nbox,scheme%grain,energy)

      CALL CPU_TIME(T1)
      CALL TIMTXT('>>> TIME USED get NF energy', T1-T0, LUPRI)

!      CALL sefmm_get_J_matrix(basis%nbas,npts,nbox,scheme%grain,Jmat)
!      CALL CPU_TIME(T1)
!      CALL TIMTXT('>>> TIME USED get NF J-matrix', T1-T0, LUPRI)

      CALL sefmm_free_potential

   END SUBROUTINE sefmm_get_NF_J_matrix

!-------------------------------------------------------------------------------

END MODULE sefmm_driver
